Fitting meta-analysis with the glmmTMB R package

Coralie Williams, Maeve McGillycuddy, David Warton, Shinichi Nakagawa

Overview

• Motivation for meta-analysis
• Model approaches
• Implementations
• Simulation study
• Illustrative example

Rapid growth of scientific research

• Doubles every ~17 years (Bornmann et al., 2021)

• Conflicting results across related studies

• Hard to keep track of the overall body of evidence

There is a lot of empirical research being published every year. A recent study found that, scientific output has been growing and doubles every 17 years. With that constant growth, it becomes very hard to keep track of the evidence and to see the bigger picture—especially when different studies report conflicting results or use different methods and measures.

Meta-analysis

Conventional meta-analysis model

1. Derive effect-size estimates and sampling variances

e.g. log-transformation or standardisation of raw data from each study

2. Fit model

Linear mixed model (LMM) with known sampling variances

\[ y_i = \mu + u_i + e_i \] \[ e_i \sim N(0, v_i)\] \[u_i \sim N(0, \tau^2) \]

For each study \(i\):

• \(y_i\) is the observed effect size.

• \(\mu\) is the overall mean effect size.

• \(u_i\) are study-level random effects.

• \(e_i\) are sampling errors with known variances \(v_i\).

Also called Normal-Normal model

Another approach

Analyse raw data using a Generalized Linear Mixed Model (GLMM) framework.

Example:

Binary comparative outcomes (e.g. treatment success/failure)

\[ X_{ij} \sim \text{Binomial}(n_{ij}, p_{ij}) \] \[ \operatorname{logit}(p_{ij}) = \mu + u_i \] \[ u_i \sim N(0, \tau^2) \]

For study \(i\) and group \(j\):

• \(X_{ij}\) is the number of events
  
• \(n_{ij}\) is the number of subjects
  
• \(p_{ij}\) is the event probability
  

Binomial–Normal model.

Two broad approaches

‘Two-stage’ approach

• Conventional meta-analysis
  
• LMM framework
• Sampling variances are known/estimated

‘One-stage’ approach

• GLMM framework
  
• Uses binary or count data from studies directly

• First stage: estimate the effect sizes and their sampling variances
• Second stage: assume normality of the empirical effect sizes
• Uses the raw data directly

Implementations

Meta-analysis packages in R

Currently > 180 packages on CRAN task view for meta-analysis

• metafor (Viechbauer, 2010): comprehensive and flexible modelling

• meta (Balduzzi et al., 2019): basic modelling and in user-friendly interface

glmmTMB R package

Increasing community > 50k downloads per month (as of 2025)

• Estimation using Template Model Builder (TMB) and Laplace approximation

• Flexible formula interface

• Range of distributional families

• Dispersion and zero-inflation modelling

➜ Currently we can fit one-stage GLMM meta-analysis models ✅

➜ But lacks functionality for two-stage models ❌

New covariance structure equalto

Random-effect term to fix the sampling covariance structure:

equalto(0 + id|g, V)

Three components:

1. id, effect size indices as a factor.

2. g, grouping variable ➜ for meta-analysis we fix this to one level.

3. V, user-supplied sampling variance–covariance matrix.

Why another implementation?

• Expands glmmTMB toolkit for meta-analysis users.

• Broader range of models with some advantages:

 ○ Complex covariance structures (Kristensen & McGillycuddy, 2025).

 e.g. reduce-rank approaches, time series, spatial structures.

 ○ Flexible dispersion formula with random effects.

 ○ Fast(er) estimation for complex models.

Simulation study

Benchmark with rma.uni function from metafor package.

Assuming random-effects model: \[ y_i = \mu + u_i + e_i, \quad u_i \sim N(0, \tau^2) \]

• Standardised mean difference (SMD)
  
• Log response ratio (lnRR)
  
• Log odds ratio (OR)
  
• Log incidence rate ratio (IRR)
  

➜ Compare estimates of the overall mean \(\hat\mu\) and within study variance \(\hat\tau^2\)

Registration of study design: https://osf.io/bvnug

Agreement between glmmTMB and metafor

Illustrative example

Association between maternal size, and number of offspring.

Dataset

Available via the metadat package (White et al., 2022): dat.lim2014

  study             species     ri  ni      yi     vi id 
1     1 Sceloporus_virgatus  0.100  21  0.1003 0.0556  1 
2     1 Sceloporus_virgatus -0.170  14 -0.1717 0.0909  2 
3     1 Sceloporus_virgatus -0.070  21 -0.0701 0.0556  3 
4     1 Sceloporus_virgatus -0.140  14 -0.1409 0.0909  4 
5     2     Marmota_marmota -0.540  74 -0.6042 0.0141  5 
6     3      Vipera_ursinii  0.487 105  0.5321 0.0098  6 

• ri correlation coefficients

• ni: sample sizes

• yi effect sizes estimates (Fisher’s z-transformed correlation coefficients)

• vi sampling variance estimates

Multiple effect size observed per study. ➜ 170 effect sizes across 120 species and 125 studies.

Phylogenetic multilevel meta-analysis

\[ y_{ijk} = \mu + u_{i} + m_{ij} + s_k + p_k + e_{ij} \]

Where:

• \(y_{ijk}\) is the observed effect size \(j\) from study \(i\) for species \(k\)

• \(u_i\) is the among-study effect and \(m_{ij}\) is within-study effect

Two species-level random effects (Cinar et al., 2022):

• \(s_k\) independent species effect

• \(p_k \sim N(0, \sigma^2_p A)\) where \(A\) is the phylogenetic correlation matrix

• \(e_{ij} \sim N(0, V)\) where \(V\) is the sampling error variance-covariance matrix.

Sampling error variance-covariance matrix \(V\):

V <- diag(dat$vi)
round(V[1:5,1:5], 3)

      [,1]  [,2]  [,3]  [,4]  [,5]
[1,] 0.056 0.000 0.000 0.000 0.000
[2,] 0.000 0.091 0.000 0.000 0.000
[3,] 0.000 0.000 0.056 0.000 0.000
[4,] 0.000 0.000 0.000 0.091 0.000
[5,] 0.000 0.000 0.000 0.000 0.014

Phylogenetic correlation matrix \(A\):

                      Hogna_helluo Centruroides_vittatus Gambusia_holbrooki
Hogna_helluo             1.0000000             0.9915966                  0
Centruroides_vittatus    0.9915966             1.0000000                  0
Gambusia_holbrooki       0.0000000             0.0000000                  1

Fit multilevel model

library(glmmTMB)
fit.tmb <- glmmTMB(yi ~ 1 + (1|study) +
                       (1 | species) + propto(0 + species.phy|g, A) +
                       equalto(0 + id|g,V),
                       data=dat,
                       REML=TRUE) 

• Use the propto covariance structure to specify phylogenetic random effect.

• Use the equalto covariance structure to specify sampling error structure.

• The within-study variability \(m_{ij}\) is modelled as the residual variance.

Model output

Fixed effects estimates from conditional model:

summary(fit.tmb)$coefficients$cond

              Estimate Std. Error   z value  Pr(>|z|)
(Intercept) -0.1563723  0.1277383 -1.224161 0.2208913

Random effect variance estimates:

vc <- VarCorr(fit.tmb)$cond
s2.m <- sigma(fit.tmb)^2 #m_ij var estimate
data.frame(
  re = c("study", "ef", "species", "species.phy"),
  var.est = c(vc$study[1], s2.m, vc$species[1], vc$g[1])
)

           re      var.est
1       study 1.387349e-01
2          ef 9.255962e-03
3     species 2.505568e-10
4 species.phy 5.721480e-02

Conclusions

• equalto supports the two-stage meta-analysis approach in glmmTMB

• Expands the meta-analysis toolkit in R

• Broadens the modelling possibilities for meta-analysis users

Next steps

• Working to include equalto on glmmTMB main branch.

• Preparing tutorial paper

Current branch: https://github.com/coraliewilliams/glmmTMB/tree/equalto_covstruc

Thank you

Supervisors:
Shinichi Nakagawa
David Warton

Maeve McGillycuddy — glmmTMB
Mollie Brooks — glmmTMB
Ben Bolker — glmmTMB
Wolfgang Viechtbauer — metafor
Yefeng Yang
Ayumi Mizuno

References

Balduzzi S, Rücker G, Schwarzer G (2019), How to perform a meta-analysis with R: a practical tutorial, Evidence-Based Mental Health; 22: 153-160.

Bornmann, L., Haunschild, R., & Mutz, R. (2021). Growth rates of modern science: A latent piecewise growth curve approach to model publication numbers from established and new literature databases. Humanities and Social Sciences Communications, 8(1), 224. https://doi.org/10.1057/s41599-021-00903-w

Brooks, M. E., Kristensen, K., van Benthem, K. J., Magnusson, A., Berg, C. W., Nielsen, A., Skaug, H. J., Mächler, M., & Bolker, B. M. (2017). glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. R Journal, 9(2), 378–400. https://doi.org/10.32614/RJ-2017-066

Cinar, O., Nakagawa, S., & Viechtbauer, W. (2022). Phylogenetic multilevel meta-analysis: A simulation study on the importance of modelling the phylogeny. Methods in Ecology and Evolution, 13(2), 383–395. https://doi.org/10.1111/2041-210X.13760

Kristensen, K., & McGillycuddy, M. (2025). Covariance structures with glmmTMB. https://cran.r-project.org/web/packages/glmmTMB/vignettes/covstruct.html

Lim, J. N., Senior, A. M., & Nakagawa, S. (2014). Heterogeneity in individual quality and reproductive trade-offs within species. Evolution, 68(8), 2306–2318. https://doi.org/10.1111/evo.12446

Viechtbauer, W. (2010). Conducting Meta-Analyses in R with the metafor Package. Journal of Statistical Software, 36(3), 1-48. https://doi.org/10.18637/jss.v036.i03

White T, Noble D, Senior A, Hamilton W, Viechtbauer W (2022). metadat: Meta-Analysis Datasets. R package version 1.2-0, https://CRAN.R-project.org/package=metadat.